En este projecto usaremos el fichero de datos “scrubbed.csv” obtenido en https://www.kaggle.com/datasets/NUFORC/ufo-sightings. Elimininaremos datos de tipo NA y cambiaremos paises vacios (i.e. pais == ““) por”Otros”. Tambien arreglamos el formato de los datos con la libreria lubridate. Usamos la libreria dplyr para eliminar columnas no necesarias.
# Leemos datos de CSV
datos <- read.csv("scrubbed.csv")
# Reemplacamos paises "" por "Otros"
datos$country <- ifelse(datos$country=="","otros",datos$country)
# Formateamos datos
datos$datetime <- as.POSIXct(datos$datetime, format = "%m/%d/%Y %H:%M")
datos$'duration (seconds)' <- as.integer(datos$'duration..seconds.')
datos$latitude <- as.numeric(datos$latitude)
# Eliminamos datos NA
datos <- na.omit(datos)
library(dplyr)
datos <- datos %>% select(-one_of('duration..seconds.', 'duration..hours.min.', 'date.posted'))
Analisando los datos por pais, lógicamente podemos ver que la mayoria de los OVNIs reportados a la National UFO Reporting Center de EEUU ocurririon dentro de EEUU. Usaremos la libreria ggplot2 para generar un diagrama representando el numero de OVNIs reportados por pais y la libreria leaflet para generar un mapa con las longitudes y latitudes.
library(leaflet)
library(ggplot2)
by_country <- aggregate(datos$country, by=list(Country = datos$country), FUN=length)
ggplot(by_country, aes(x=by_country$Country, y=by_country$x, fill=by_country$Country)) +
geom_bar(stat="identity", width=1) +
xlab('Paises') +
ylab('Avistamientos') +
guides(fill = guide_legend(title = "Country" ))
leaflet(datos) %>%
addTiles() %>%
addProviderTiles("CartoDB.Positron") %>%
addMarkers(~longitude, ~latitude,
popup = ~state, label = ~city,
clusterOptions = markerClusterOptions())
Por eso, hemos decidido hacer nuestros analisis sobre datos solamente de EEUU. Tomamos una prueba de tamaño 1000.
datos_us <- datos[which(datos$country == 'us'),]
datos_us <- datos_us %>% select(-country)
datos_us <- datos_us[sample(nrow(datos_us), 1000), ]
De la prueba de 1000 datos, se ha hecho el análisis solamente de las variables cuantitativas, las cuales son: (“datetime, latitude, longitude, duration(seconds)”).
Con el comando names(object…) se puede ver todas las columnas del fichero de datos y las referenciadas anteriormente.
names(datos_us)
## [1] "datetime" "city" "state"
## [4] "shape" "comments" "latitude"
## [7] "longitude" "duration (seconds)"
Con summary(object...) podemos ver la media, mediana y
percentiles de varios datos, para calcular la moda usamos la funcion
mode.
mode <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
summary(datos_us$datetime)
## Min. 1st Qu. Median
## "1947-04-15 14:00:00" "2001-11-03 17:36:15" "2007-07-24 05:37:30"
## Mean 3rd Qu. Max.
## "2004-10-25 21:59:05" "2011-09-09 05:12:30" "2014-05-06 22:30:00"
print(paste("Moda: ", mode(datos_us$datetime)))
## [1] "Moda: 2012-07-04 22:30:00"
summary(datos_us$latitude)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 19.73 34.05 38.74 38.22 41.85 61.22
print(paste("Moda: ", mode(datos_us$latitude)))
## [1] "Moda: 33.4483333"
summary(datos_us$longitude)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -158.19 -114.91 -91.64 -96.28 -81.03 -68.41
print(paste("Moda: ", mode(datos_us$longitude)))
## [1] "Moda: -112.0733333"
summary(datos_us$'duration (seconds)')
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 30 180 1086 600 109800
print(paste("Moda: ", mode(datos_us$'duration (seconds)')))
## [1] "Moda: 300"
No se ha considerado necesario calcular la varianza y desviación poblacional
| Medida | Formula |
|---|---|
| Cuasivarianza | var(x) |
| Cuasidesviación típica | sd(x) |
| Coeficiente de variación | coef_var(x) |
coef_var <- function(x) {
sd(x) / mean(x)
}
varianza <- var(datos_us$'duration (seconds)')
print(paste("Cuasivarianza: ", varianza))
## [1] "Cuasivarianza: 35067587.5226226"
sd <- sd(datos_us$'duration (seconds)')
print(paste("Cuasidesviación típica: ", sd))
## [1] "Cuasidesviación típica: 5921.78921632834"
cv <- coef_var(datos_us$'duration (seconds)')
print(paste("Coeficiente de variación: ", cv))
## [1] "Coeficiente de variación: 5.4511881438682"
varianza <- var(datos_us$latitude)
print(paste("Cuasivarianza: ", varianza))
## [1] "Cuasivarianza: 29.2851452840036"
sd <- sd(datos_us$latitude)
print(paste("Cuasidesviación típica: ", sd))
## [1] "Cuasidesviación típica: 5.41157512042506"
cv <- coef_var(datos_us$latitude)
print(paste("Coeficiente de variación: ", cv))
## [1] "Coeficiente de variación: 0.141580908541579"
varianza <- var(datos_us$longitude)
print(paste("Cuasivarianza: ", varianza))
## [1] "Cuasivarianza: 326.357120221342"
sd <- sd(datos_us$longitude)
print(paste("Cuasidesviación típica: ", sd))
## [1] "Cuasidesviación típica: 18.0653569082192"
cv <- coef_var(datos_us$longitude)
print(paste("Coeficiente de variación: ", cv))
## [1] "Coeficiente de variación: -0.187626786221631"
Como no se pueden dividir fechas, no calculamos el Coeficiente de variación para datetime.
varianza <- var(datos_us$datetime)
print(paste("Cuasivarianza: ", varianza))
## [1] "Cuasivarianza: 108644846179268224"
sd <- sd(datos_us$datetime)
print(paste("Cuasidesviación típica: ", sd))
## [1] "Cuasidesviación típica: 329613176.586235"
Al calcular la asimetría y curtosis de la latitud determinamos que: Obtenemos una asimetría negativa, por tanto la mayoría de datos se encuentran a la derecha de la media. Obtenemos una curtosis mayor que cero, por tanto, hay una mayor concentración de datos alrededor de la media
library(moments)
x <- (datos_us$latitude)
print(kurtosis(x))
## [1] 3.553226
print(skewness(x))
## [1] -0.04729688
print(mean(x))
## [1] 38.22249
hist(x)
### longitude Al calcular la asimetría y curtosis de la latitud
determinamos que: Obtenemos una asimetría negativa, por tanto la
mayoría de datos se encuentran a la derecha de la media. Obtenemos
una curtosis mayor que cero, por tanto, hay una mayor concentración de
datos alrededor de la media
library(moments)
x <- (datos_us$longitude)
print(kurtosis(x))
## [1] 2.315487
print(skewness(x))
## [1] -0.4947675
print(mean(x))
## [1] -96.28346
hist(x)
### datetime Al calcular la asimetría y curtosis de la latitud
determinamos que: Obtenemos una asimetría negativa, por tanto la
mayoría de datos se encuentran a la derecha de la media. Obtenemos
una curtosis mayor que cero, por tanto, hay una mayor concentración de
datos alrededor de la media
library(moments)
x <- (datos_us$datetime)
print(kurtosis(x))
## [1] 10.29702
print(skewness(x))
## [1] -2.502833
hist(x, 50)
### duration (seconds) Al calcular la asimetría y curtosis de la latitud
determinamos que: Obtenemos una asimetría positiva, por tanto la
mayoría de datos se encuentran a la izquierda de la media.
Obtenemos una curtosis mayor que cero, por tanto, hay una mayor
concentración de datos alrededor de la media
library(moments)
x <- (datos_us$'duration (seconds)')
print(kurtosis(x))
## [1] 264.1936
print(skewness(x))
## [1] 15.34129
print(mean(x))
## [1] 1086.33
Asociados a la Tabla de Frecuencias para algunos datos (excepto histogramas, usados anteriormente para ver mejor las medidas de forma) ### latitud
breaks <- seq(0, 90, by=10)
lat.cumfreq <- c(0,cumsum(table(cut(datos_us$latitude, breaks, right = FALSE))))
plot(breaks, lat.cumfreq,
main="Latitud de los ovnis",
xlab = "Latitud",
ylab = "Latitudes acumuladas")
lines(breaks, lat.cumfreq)
### longitude
breaks <- seq(-200, 100, by=10)
longitude.cumfreq <- c(0,cumsum(table(cut(datos_us$longitude, breaks, right = FALSE))))
plot(breaks, longitude.cumfreq,
main="Longitud de los ovnis",
xlab = "Longitud",
ylab = "Longitudes acumuladas")
lines(breaks, longitude.cumfreq)
Si quieres tratar de avistar un UFO la siguiente información que te proporcionamos te será de suma utilidad, como el pais, estado y ciudad donde debes ir para realizar un avistamiento, además de la forma que debes buscar y las horas a las que debes buscarlo, además de otros datos curiosos como la duración y los años que fueron más probables de avistar.
print("Probabilidad de avistar un ovni por país:")
## [1] "Probabilidad de avistar un ovni por país:"
dat <- read.csv("scrubbed.csv")
dat_country<-aggregate(dat$country, by=list(Country = dat$country), FUN=length)
dat$country <- ifelse(dat$country=="","otros",dat$country)
numAvistamientos<-sum(dat_country$x)
proobau<-dat_country[which(dat_country$Country == 'au'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Australia: ",probbauval,"%"))
## [1] "Australia: 0.669720659264054 %"
proobau<-dat_country[which(dat_country$Country == 'ca'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Canadá: ",probbauval,"%"))
## [1] "Canadá: 3.73450181745755 %"
proobau<-dat_country[which(dat_country$Country == 'de'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Alemania: ",probbauval,"%"))
## [1] "Alemania: 0.130707563611014 %"
proobau<-dat_country[which(dat_country$Country == 'gb'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Gran Bretaña: ",probbauval,"%"))
## [1] "Gran Bretaña: 2.37140865408555 %"
proobau<-dat_country[which(dat_country$Country == 'us'),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("Estados Unidos: ",probbauval,"%"))
## [1] "Estados Unidos: 81.056117113977 %"
proobau<-dat_country[which(dat_country$Country == ''),]
probbauval<-(proobau$x/numAvistamientos)*100
print(paste("otros países: ",probbauval,"%"))
## [1] "otros países: 12.0375441916048 %"
maxcountry<-sum(by_country$x)
ggplot(by_country, aes(x=Country, y=(x/maxcountry)*100, fill=Country)) +
geom_bar(stat="identity", width=1) +
xlab('Paises') +
ylab('Probabilidad') +
guides(fill = guide_legend(title = "Country" ))
Debido a que es en Estados Unidos donde hay mayor probabilidad de avistar un UFO vamos a calcular los estados donde es mas probable lograrlo
datos_us$state <- ifelse(datos_us$state=="","otros",datos_us$state)
by_state <- aggregate(datos_us$state, by=list(State = datos_us$state), FUN=length)
avistTot<-sum(by_state$x)
mostseekstate<-by_state[which(by_state$x>16),]
maxstate<-by_state[which(by_state$x == max(by_state$x)),]
print(paste("El estado con mayor probabilidad de avistar un ovni es",maxstate$State,"con una probabilidad de un ",(maxstate$x/avistTot)*100,"%"))
## [1] "El estado con mayor probabilidad de avistar un ovni es ca con una probabilidad de un 13.8 %"
print(paste("A continuación se muestra una gráfica con los estados donde mayor probabilidad hay de avistar un ovni"))
## [1] "A continuación se muestra una gráfica con los estados donde mayor probabilidad hay de avistar un ovni"
ggplot(mostseekstate, aes(x=mostseekstate$State, y=(mostseekstate$x/avistTot)*100, fill=mostseekstate$State)) +
geom_bar(stat="identity", width=1) +
xlab('') +
ylab('Probabilidad') +
theme(axis.text.x=element_blank())+
guides(fill = guide_legend(title = "Estados" ))
Ya que hemos calculado los estados haremos lo mismo con las ciudades
by_city <- aggregate(datos_us$city, by=list(City = datos_us$city), FUN=length)
avistTot<-sum(by_city$x)
mostseek<-by_city[which(by_city$x>4),]
maxcity<-by_city[which(by_city$x == max(by_city$x)),]
print(paste("La ciudad con mayor probabilidad de avistar un ovni es",maxcity$City,"con una probabilidad de un ",(maxcity$x/avistTot)*100,"%"))
## [1] "La ciudad con mayor probabilidad de avistar un ovni es phoenix con una probabilidad de un 1.2 %"
print(paste("A continuación se muestra una gráfica con las ciudades donde mayor probabilidad hay de avistar un ovni"))
## [1] "A continuación se muestra una gráfica con las ciudades donde mayor probabilidad hay de avistar un ovni"
ggplot(mostseek, aes(x=City, y=(x/avistTot)*100, fill=City)) +
geom_bar(stat="identity", width=1) +
xlab('') +
ylab('Probabilidad') +
theme(axis.text.x=element_blank())+
guides(fill = guide_legend(title = "Ciudades" ))
En cuanto a la forma que tienes que buscar en el cielo para lograr un avistamiento exitoso
by_shape <- aggregate(datos_us$shape, by=list(shape = datos_us$shape), FUN=length)
by_shape$shape <- ifelse(by_shape$shape=="","otros",by_shape$shape)
max_shape<-sum(by_shape$x)
maxshape<-by_shape[which(by_shape$x == max(by_shape$x)),]
minshape<-by_shape[which(by_shape$x == min(by_shape$x)),]
print(paste("El ovni mas probabilidad de ser avistado son aquellos ovnis con forma",maxshape$shape,"con una probabilidad del",(maxshape$x/avistTot)*100,"%"))
## [1] "El ovni mas probabilidad de ser avistado son aquellos ovnis con forma light con una probabilidad del 21.2 %"
print(paste("Mientras que los más raros con mayor dificultad de avistamiento son aquellos con forma"))
## [1] "Mientras que los más raros con mayor dificultad de avistamiento son aquellos con forma"
print (paste(minshape$shape))
## [1] "cross"
print(paste("con una probabilidad de avistamiento del",(minshape$x[1]/avistTot)*100,"%"))
## [1] "con una probabilidad de avistamiento del 0.2 %"
ggplot(by_shape, aes(x=by_shape$shape,y=(by_shape$x/max_shape)*100,fill=by_shape$shape)) +
geom_bar(stat="identity", width=1) +
xlab('') +
ylab('Avistamientos') +
theme(axis.text.x=element_blank())+
guides(fill = guide_legend(title = "Formas" ))
La hora también es importante, hay que saber cuando buscar y cuando
descansar
by_horas<- aggregate(datos_us$datetime, by=list(Hour = format(datos_us$datetime,"%H")), FUN=length)
maxhora<-by_horas[which(by_horas$x == max(by_horas$x)),]
maxhour<-sum(by_horas$x)
print(paste("La hora en la que es más probable avistar un ovni son las",maxhora$Hour,':00',"con una probabilidad del",(maxhora$x/avistTot)*100,"%"))
## [1] "La hora en la que es más probable avistar un ovni son las 22 :00 con una probabilidad del 14.4 %"
plot(by_horas$Hour,(by_horas$x/maxhour)*100,type="s",col="dark red", xlab = "Horas del día", ylab = "Probabilidad", main = "Probabilidad de avistamiento por hora")
También es interesante ver cuando hubo una mayor probabilidad de avistar un ovni
anios<- aggregate(datos_us$datetime, by=list(Date = format(datos_us$datetime,"%Y")), FUN=length)
maxanio<-anios[which(anios$x == max(anios$x)),]
maxyear<-sum(anios$x)
print(paste("El año con mayor probabilidad en caso de avistamiento fué",maxanio$Date,"con una probabilidad del",(maxanio$x/avistTot)*100,"%"))
## [1] "El año con mayor probabilidad en caso de avistamiento fué 2012 con una probabilidad del 9.9 %"
plot(anios$Date,(anios$x/maxyear)*100,type="l",col="dark blue", xlab = "Linea temporal", ylab = "Probabilidad", main = "Probabilidad de avistamientos por año")
Observamos la probabilidad de lograr un avistamiento de larga
duración
durability<- aggregate(datos_us$`duration (seconds)`, by=list(duracion = datos_us$`duration (seconds)`), FUN=length)
maxduracion<-sum(durability$x)
maxdur<-durability[which(durability$x>1),]
durationProb<-durability[which(durability$x<=10),]
sumaprobduracion<-sum(durationProb$x)
print(paste("La probabilidad de ver un UFO durante menos de 10 segundos es del",(sumaprobduracion/maxduracion)*100,"%"))
## [1] "La probabilidad de ver un UFO durante menos de 10 segundos es del 12.3 %"
durationProb<-durability[which(durability$x<=60),]
sumaprobduracion<-sum(durationProb$x)
print(paste("La probabilidad de ver un UFO durante menos de un minuto es del",(sumaprobduracion/maxduracion)*100,"%"))
## [1] "La probabilidad de ver un UFO durante menos de un minuto es del 60.9 %"
durationProb<-durability[which(durability$x>60),]
sumaprobduracion<-sum(durationProb$x)
print(paste("La probabilidad de ver un UFO durante mas de un minuto es del",(sumaprobduracion/maxduracion)*100,"%"))
## [1] "La probabilidad de ver un UFO durante mas de un minuto es del 39.1 %"
plot(maxdur$duracion,(maxdur$x/sumaprobduracion)*100,type="l",col="blue", xlab = "Segundos", ylab = "Probabilidad", main = "Probabilidad de avistamientos por más de un segundo")
A continuación realizaremos algunas operaciones con sucesos de
probabilidad
prob1 <- datos_us[which(datos_us$city == 'seattle'&format(datos_us$datetime,"%H")>12),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en seattle por la tarde",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en seattle por la tarde 0.5 %"
prob1 <- datos_us[which(datos_us$city == 'seattle'|format(datos_us$datetime,"%H")==21),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en seattle o verlo a las 21:00 es del",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en seattle o verlo a las 21:00 es del 13.7 %"
prob1 <- datos_us[which(datos_us$city != 'seattle'&format(datos_us$datetime,"%H")!=21),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en una ciudad que no sea seattle a una hora distinta a las 21:00 es del",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en una ciudad que no sea seattle a una hora distinta a las 21:00 es del 86.3 %"
prob1 <- datos_us[which(datos_us$city != 'seattle'&format(datos_us$datetime,"%H")==21),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en una ciudad que no sea seattle a las 21:00 es del",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en una ciudad que no sea seattle a las 21:00 es del 13.2 %"
prob1 <- datos_us[which(datos_us$city != 'seattle'&format(datos_us$datetime,"%H")==21&(datos_us$shape == 'light'|datos_us$shape=='unknown')),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en una ciudad distinta de seattle a las 21:00 con una forma que sea lumínica o desconocida",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en una ciudad distinta de seattle a las 21:00 con una forma que sea lumínica o desconocida 4 %"
prob1 <- datos_us[which(datos_us$city == 'seattle'&format(datos_us$datetime,"%H")!=21&(datos_us$shape != 'light'&datos_us$shape!='unknown')&datos_us$`duration (seconds)`>10),]
prob2<-aggregate(prob1$city,by=list(city=prob1$city),FUN=length)
print(paste("La probabilidad de ver un UFO en seattle a una hora distinta de las 21:00 con una forma que sea distinta de lumínica o desconocida pero con una duracion de mas de 10 segundos",(sum(prob2$x)/avistTot)*100,"%"))
## [1] "La probabilidad de ver un UFO en seattle a una hora distinta de las 21:00 con una forma que sea distinta de lumínica o desconocida pero con una duracion de mas de 10 segundos 0.3 %"
Respondemos algunas preguntas de probabilidad condicional
maxdat<-count(datos)
prob1 <- datos_us[which(format(datos_us$datetime,"%H")==15),]
prob2<-count(prob1)
print(paste("Dado que son las 15:00 que probabilidad hay de ver un UFO?",(prob2/avistTot)*100,"%"))
## [1] "Dado que son las 15:00 que probabilidad hay de ver un UFO? 1.6 %"
prob1 <- datos_us[which(datos_us$city == 'richmond'),]
prob2<-count(prob1)
print(paste("Si vivo en Richmond que probabilidad tengo de ver un UFO?",(prob2/avistTot)*100,"%"))
## [1] "Si vivo en Richmond que probabilidad tengo de ver un UFO? 0.3 %"
prob1 <- datos_us[which(datos_us$shape=='changing'),]
prob2<-count(prob1)
print(paste("En caso de ver un ovni que probabilidad hay de que sea cambiante?",(prob2/avistTot)*100,"%"))
## [1] "En caso de ver un ovni que probabilidad hay de que sea cambiante? 2.5 %"
prob1 <- datos[which(datos$country != 'us'),]
prob2<-count(prob1)
print(paste("Si vivo fuera de Estados Unidos que probabilidad hay de avistar un ovni?",(prob2/maxdat)*100,"%"))
## [1] "Si vivo fuera de Estados Unidos que probabilidad hay de avistar un ovni? 18.9425722359854 %"
Como los datos, excepto los extremos, están aproximadamente sobre la linea, podemos decir que exluyendo los extremos, los datos siguen una distribucion normal.
dia <-as.numeric(format(datos_us$datetime,'%d'))
qqnorm(dia, pch = 1, frame = FALSE)
qqline(dia, col = "steelblue", lwd = 2)
También tenemos la media, moda, varianza y desviación típica de estos
datos.
mean(dia)
## [1] 14.984
mode(dia)
## [1] 15
var(dia)
## [1] 79.80155
sd(dia)
## [1] 8.933171
Funcion de distribucion acumulada respecto a los días de los mes mas comunes de todos los años en el que se han visto ovnis en EEUU.
f1 <- pnorm(dia,mean(dia),sd(dia))
plot(dia,f1,main = "Funcion de distribucion acumulada")
Calculamos la probabilidad de que se avisten ovnis durante los primeros
10 dias de cada mes
P(X≤10)
pnorm(10,mean(dia),sd(dia))
## [1] 0.2884493
Calculamos la probabilidad de que se avisten ovnis durante los ultimos 5 dias de cada mes
P(X>25)
1-pnorm(25,mean(dia),sd(dia))
## [1] 0.1310983
Calculamos la probabilidad de que se avisten ovnis entre los dias 10 y 20 de cada mes.
P(10≤X≤20)
pnorm(20,mean(dia),sd(dia)) - pnorm(10,mean(dia),sd(dia))
## [1] 0.4243233
Podemos representar este último gráficamente:
regionX=seq(10,20,0.01)
xP <- c(10,regionX,20)
yP <- c(0,dnorm(regionX,15,12),0)
curve(dnorm(x,15,12),xlim=c(0,30),yaxs="i",ylim=c(0,0.035),ylab="f(x)",
main='Densidad P(10<X<20)')
polygon(xP,yP,col="orange1")
box()
Podemos averiguar cuantos ovnis fueron divisados en 2012.
year <-as.numeric(format(datos_us$datetime,'%y'))
landa <- sum(year == 12)
plot(dpois(0:200, landa), type = "h", lwd = 2,
main = "Función de masa de probabilidad",
ylab = "P(X = x)", xlab = "Número de ovnis avistados en 2012")
Una vez calculada la tasa de llegada de los ovnis en 2012 con la distribucion de Poisson, podemos dibujar la gráfica de la función de densidad exponencial.
# Rejilla del eje X
x <- seq(0, 0.1, 0.01)
dexp(x,landa)
## [1] 99.000000000 36.786092411 13.668854494 5.079027723 1.887248315
## [6] 0.701257484 0.260570935 0.096822086 0.035976830 0.013368152
## [11] 0.004967294
# lambda
plot(x, dexp(x, landa), type = "l",
ylab = "f(x)", lwd = 2, col = "red")
De la misma forma tenemos la funcion de distribución de exponencial
acumulada
# Rejilla de valores del eje X
x <- seq(0,0.1, 0.01)
pexp(x,landa)
## [1] 0.0000000 0.6284233 0.8619308 0.9486967 0.9809369 0.9929166 0.9973680
## [8] 0.9990220 0.9996366 0.9998650 0.9999498
# lambda
plot(x, pexp(x, landa), type = "l",
ylab = "F(x)", lwd = 2, col = "orange")
La probabilidad de que se divisen menos de 100 ovnis en 2012 P(X<=100) es:
ppois(100,landa)
## [1] 0.5663565
Probabilidad de que se divisen entre 80 y 95 ovnis en 2012 P(80< X< 95)
A modo ilustrativo, podemos representarlo en una gráfica
ppois(95,landa) - ppois(80,landa)
## [1] 0.3396651
pois_sum <- function(lambda, lb, ub, col = 4, lwd = 1, ...) {
x <- 0:(lambda + lambda * 2)
if (missing(lb)) {
lb <- min(x)
}
if (missing(ub)) {
ub <- max(x)
}
plot(dpois(x, lambda = lambda), type = "h", lwd = lwd, ...)
if(lb == min(x) & ub == max(x)) {
color <- col
} else {
color <- rep(1, length(x))
color[(lb + 1):ub ] <- col
}
lines(dpois(x, lambda = lambda), type = "h",
col = color, lwd = lwd, ...)
}
pois_sum(landa,80, 95, lwd = 2,
col = 2, ylab = "P(X = x)", xlab = "Ovnis avistados en 2012")
Gráfico de la función de distribución exponencial
plot(ppois(40:100, landa), type = "s", lwd = 2,
main = "Función de distribución exponencial",
xlab = "Número de eventos", ylab = "F(x)")
También podemos obtener los cuantiles correspondientes de la distribucion exponencial
plot(qpois(seq(0,1,0.001), landa),
main = "Función cuantil",
ylab = "Q(p)", xlab = "Cuantiles",
type = "s", col = 3, xaxt = "n")
axis(1, labels = seq(0, 1, 0.1), at = 0:10 * 100)
Para los datos, vamos a calcular la probabilidad de se vean ovnis el dia 24 de cada mes.
Con lo que si se ven ovnis los días 24, son aciertos, y al contrario son fallos.
dia <-as.numeric(format(datos_us$datetime,'%d'))
tam<- 1000
x <- (dia == 24)
y <- sum(x)
probEnsayo <- (y/tam)
print(probEnsayo)
## [1] 0.03
Funcion de probabilidad binomial
plot(dbinom(1:100, tam, probEnsayo), type = "h", lwd = 2,
main = "Función de probabilidad binomial",
ylab = "P(X = x)", xlab = "Número de éxitos")
Probabilidad de encontrar a lo largo de 69 años (2014 - 1945) mas de 35
ovnis el día 24 de cada mes:
P(X > 35)
1 - pbinom(35, tam, probEnsayo)
## [1] 0.1539237
Probabilidad de encontrar a lo largo de 69 años (2014 - 1945) menos de 25 ovnis el día 24 de cada mes:
P(X <= 25)
pbinom(25, tam, probEnsayo)
## [1] 0.2044652
Ahora graficmente la probabilidadd de encontrar entre 25 y 35 ovnis.
P(25 <= X <= 35)
binom_sum <- function(size, prob, lb, ub, col = 4, lwd = 1, ...) {
x <- 0:size
if (missing(lb)) {
lb <- min(x)
}
if (missing(ub)) {
ub <- max(x)
}
plot(dbinom(x, size = size, prob = prob), type = "h", lwd = lwd, ...)
if(lb == min(x) & ub == max(x)) {
color <- col
} else {
color <- rep(1, length(x))
color[(lb + 1):ub ] <- col
}
lines(dbinom(x, size = size, prob = prob), type = "h",
col = color, lwd = lwd, ...)
}
binom_sum(0 :100, probEnsayo, lb = 25, ub = 35, lwd = 2,
ylab = "P(X = x)", xlab = "Número de éxitos")
Funcion de distribucion binomial
plot(pbinom(10:50, tam, probEnsayo), type = "s", lwd = 2,
main = "Función de distribución binomial",
xlab = "Número de éxitos", ylab = "F(x)")
Hemos usado la distribucion binomial para el teorema central del límite,
Realizamos el histograma de las 6 últimas muestras y de la media muestral,
tamMuestra=1000
numMuestras=25
#hemo
#Se deja ese hueco vacío aunque salga warning
mat=matrix(, nrow = numMuestras, ncol = tamMuestra)
mediax=vector()
# Genera numMuestras muestras aleatorias
# Para cada muestra, almacena la media en el vector anterior
for (i in 1:numMuestras){
x=rbinom(tamMuestra,200,probEnsayo)
mat[i,]=x
mediax[i]=mean(x)
}
# Muestra la media de las 6 últimas muestras generadas
print(c(mediax[numMuestras],mediax[numMuestras-1],mediax[numMuestras-2],
mediax[numMuestras-3],mediax[numMuestras-4],mediax[numMuestras-5]))
## [1] 6.038 6.027 5.995 5.961 5.875 5.935
# Muestra el histograma de las 6 últimas muestras y de la media muestral
par(mfrow=c(3, 3))
hist(mat[numMuestras,],probability=TRUE)
hist(mat[numMuestras-1,],probability=TRUE)
hist(mat[numMuestras-2,],probability=TRUE)
hist(mat[numMuestras-3,],probability=TRUE)
hist(mat[numMuestras-4,],probability=TRUE)
hist(mat[numMuestras-5,],probability=TRUE)
hist(mediax,probability=TRUE)
Vamos a construir un intervalo de confianza del 95% para la duracion media de los avistamientos en EEUU.
#filtramos por duracion menor o igual a 7200, para descartar valores muy altos que son inútiles en nuestro análisis
retval <- subset(datos_us, datos_us$`duration (seconds)` < 7201, select = c(state, `duration (seconds)`))
intervalo <- t.test(retval$'duration (seconds)', conf.level = 0.95)#sacamos el intervalo de confianza
print(intervalo$conf.int)#intervalo de confianza
## [1] 490.9120 624.0471
## attr(,"conf.level")
## [1] 0.95
print(mean(retval$'duration (seconds)'))#media
## [1] 557.4796
print(paste("El resultado indica que la media de la variable duracion (seconds) es de ", mean(retval$duration..seconds.), " el cual se encuentra con una confianza del 95% en el intervalo ", intervalo$conf.int[1], " ", intervalo$conf.int[2]))
## [1] "El resultado indica que la media de la variable duracion (seconds) es de NA el cual se encuentra con una confianza del 95% en el intervalo 490.912046871598 624.047136801871"
car::qqPlot(retval$'duration (seconds)', pch=19,
main='QQplot para la duración de avistamientos',
xlab='Cuantiles teóricos',
ylab='Cuantiles muestrales')
## [1] 99 140
hist(retval$'duration (seconds)', freq=TRUE,
main='Histograma para la duración de avistamientos',
xlab='Duracion (s)',
ylab='Frecuencia')
#Podemos observar que la mayoría de avistamientos duran muy poco tiempo.
Dado que la mayoría de los datos son de EEUU, queremos saber si los OVNIs realmente prefieren estar ahí. Para esto, haremos un contraste de la duración del avistamiento entre EEUU y todos los demás. Podemos asumir que como los datos son del National UFO Reporting Center (NUFORC) basado EEUU, datos de otros países no son datos poblacionales.
Basaremos el estudio en las siguentes hipótesis:
| Hipótesis | Contraste |
|---|---|
| H0 (Nula) | µ(EEUU) = µ(Otros) |
| H1 (Alternativa) | µ(EEUU) != µ(Otros) |
Donde µ(EEUU) y µ(Otros) son la duración media de un avistamiento en EEUU y en otros países respectivamente. Tomaremos el nivel de significancia α = 0.1
alpha <- 0.1
datos_us <- datos[which(datos$country == 'us'),]
datos_us <- datos_us %>% select(-country)
datos_otros <- datos[which(datos$country != 'us'),]
datos_otros <- datos_otros %>% select(-country)
# Eliminamos valores extremos
duracion_us <- datos_us$`duration (seconds)`[!(datos_us$`duration (seconds)` %in% boxplot.stats(datos_us$`duration (seconds)`)$out)]
duracion_otros <- datos_otros$`duration (seconds)`[!(datos_otros$`duration (seconds)` %in% boxplot.stats(datos_otros$`duration (seconds)`)$out)]
Como tenemos un número de datos muy grande, usamos el teorema central del límite (verificando que es distribución normal) para calcular µ(EEUU). Dado que las duraciones son muy variables, asumimos que las desviaciones desconocidas son distintas.
muestra_duracion_us <- replicate(100, sample(duracion_us, size = 1000, replace = TRUE) %>% mean())
car::qqPlot(muestra_duracion_us)
## [1] 70 97
var_us <- var(muestra_duracion_us)
x_us <- mean(muestra_duracion_us)
Hacemos lo mismo para las duraciones de otros países
muestra_duracion_otros <- replicate(100, sample(duracion_otros, size = 1000, replace = TRUE) %>% mean())
car::qqPlot(muestra_duracion_otros)
## [1] 96 92
var_otros <- var(muestra_duracion_otros)
x_otros <- mean(muestra_duracion_otros)
Calculamos la region de rechazo
f <- ((var_us / 100 + var_otros / 100) ^ 2) / ((var_us/100)^2/99 + (var_otros/100)^2/99)
f <- round(f)
cuantil <- qt(c(.025, .975), df=f)
print(paste("La region de rechazo es T < ", cuantil[1], " o T > ", cuantil[2]))
## [1] "La region de rechazo es T < -1.97201747783631 o T > 1.97201747783631"
estadistico <- abs(x_us - x_otros) / sqrt(var_us / 100 + var_otros / 100)
print(estadistico)
## [1] 5.945632
Como nuestro estadistico está en la region de rechazo, rechazamos la hipótesis nula. Eso significa que efectivamente, la duración de un avistamiento es mayor en EEUU.
Podemos verificar esto usando el test t de student ya que nuestra n > 30
t.test(x=muestra_duracion_us, y=muestra_duracion_otros, alternative = "two.sided", conf.level = 0.9)
##
## Welch Two Sample t-test
##
## data: muestra_duracion_us and muestra_duracion_otros
## t = 5.9456, df = 197.71, p-value = 1.228e-08
## alternative hypothesis: true difference in means is not equal to 0
## 90 percent confidence interval:
## 6.158525 10.899955
## sample estimates:
## mean of x mean of y
## 264.8760 256.3468